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Abstract 

In this paper, we study the inverse electromagnetic medium scattering problem of estimating the 
support and shape of medium scatterers from scattered electric or magnetic near-field data. We shall 
develop a novel direct sampling method based on an analysis of electromagnetic scattering and the 
behavior of the fundamental solution. The method is applicable even with one incident field and 
needs only to compute inner products of the measured scattered field with the fundamental solutions 
located at sampling points. Hence it is strictly direct, computationally very efficient, and highly 
tolerant to the presence of noise in the data. Two- and three-dimensional numerical experiments 
indicate that it can provide reliable support estimates of one single and multiple scatterers in case of 
both exact and highly noisy data. 

Key Words: inverse medium scattering, direct sampling method, scattering analysis, electromag- 
netic wave propagation 

1 Introduction 

Inverse electromagnetic scattering represents an important noninvasive imaging technology for interrogat- 
ing material properties, and it arises in many practical applications such as biomedical diagnosis [ , ], 
nondestructive testing [9], and geophysical exploration [12]. In this work we are concerned with the inverse 
medium problem of determining the electrical/magnetic properties of unknown inhomogeneous objects 
embedded in a homogeneous background from noisy measurements of the scattered electric/magnetic 
field corresponding to one or several incident fields impinged on the objects. Mathematically, the medium 
scattering problem is described by the time-harmonic Maxwell system: 

iujcE X H = in R^, 
-luoiiH X E = in R^, 

where the vectorial fields H and E denote the magnetic and electric fields, respectively. Here the con- 
stant uj is the angular frequency, the functions e and /i refer to the electrical permittivity and magnetic 
permeability, respectively. 

Let the domain C R^ (d = 2, 3) be the space occupied by the inhomogeneous medium objects within 
the homogeneous background R^. We are interested in either electrical or magnetic inhomogeneities, but 
shall focus our discussions on the case of electrical inhomogeneities since the case of magnetic inho- 
mogeneities follows analogously. So we shall assume /i = /io, i-e., the magnetic permeability of the 
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background. Then by taking the curl of the second equation of the system and ehminating H in the first 
equation, we obtain the following vector Helmholtz equation for the electric field E 

V x{V xE)- k^n^{x) E = (1) 

where k is the wavenumber, with k'^ = cj^eo/io? and n = ^^^J^ is the refractive index function, i.e., the 
ratio of the wave velocity in the homogeneous background medium to that in the inhomogeneities. The 
refractive index n completely characterizes the electrical inhomogeneities, and the support is filled by the 
inhomogeneous media, i.e., supp(n^ — 1) = Q. Further, we assume that the medium scattering is excited 
by an incident plane wave E'^: 

E' =pe^^^-^, 

where d G S^~^ and p G S^~^ are the incident and polarization directions, respectively. Since the incident 
field E'^ is solenoidal, i.e., V-£^* = 0, the polarization p should be chosen such that it is perpendicular to the 
incident direction d. Then the incident field E^ satisfies the Maxwell system (1) in the entire homogeneous 
background. The forward scattering problem is to find the total electric field E given the refractive index 
n^, under the following Silver- Miiller radiation condition for the scattered field E^ = E — E^: 

lim \x\^ {V xE' xx-\kE')={) 

\x\^oo 

uniformly for all directions x = x/\x\ G S^~^. 

The inverse problem of our interest is the inverse medium scattering problem, which is to reconstruct 
the inhomogeneous media from the scattered electric field E^ corresponding to one (or several) incident 
field E\ measured over a certain closed curve/surface P. In view of the practical significance of the 
inverse problem, there has been considerable interest in designing efficient and stable inversion techniques. 
However, this is very challenging because of a number of complicating factors: strong nonlinearity of 
the map from the refractive index to the scattered field, severe ill-posedness of the inverse problem, 
complexity of the forward model, and the limited available and noisy data. Nonetheless, a number 
of imaging algorithms have been developed in the literature, which can be roughly divided into two 
categories: direct and indirect methods. The former aims at detecting the scatterer support and shape, 
and includes linear sampling method (LSM) [10, ], multiple signal classification (MUSIC) [ ], and 
asymptotic analysis [ ]. In contrast, the latter provides a distributed estimate of the refractive index by 
applying regularization techniques. We refer interested readers to the adjoint-based method [15, 30, 22], 
the recursive linearization (with continuation in frequency) [ , ], the Gauss- Newton method [18, 13, 19], 
the contrasted source inversion [ ], subspace regularization [ ] and level set method [16] for an incomplete 
list. Generally, the estimates by the method of the latter category can provide more details of the 
inclusions/inhomogeneities, but at the expense of much increased computational efforts, especially when 
the forward model is the full Maxwell system. 

In this work, we shall develop a novel direct method, called the direct sampling method, for stably 
and accurately detecting the support of the scatterers. These sampling-type methods are direct in the 
sense that no optimizations or solutions of linear systems are involved, so they are computationally very 
cheap. They have attracted some attentions in very recent years; see [26] and [21] for inverse acoustic 
scattering problems using far-field data and near-field data, respectively. A different derivation of the 
method due to Potthast [26] was given in [_ ], where the performance of the methods using near-field 
and far-field data was compared in detail and the effectiveness of these methods was also investigated for 
several important scattering scenarios: obstacles, inhomogeneous media, cracks and their combinations. 

The major goal of this work is to extend the direct sampling method developed in [ ] for the acoustic 
scattering to the electromagnetic scattering. Due to the much increased complexity of the Maxwell's 
system relative to its scalar counterpart, the Helmholtz equation, the extension is nontrivial and requires 
several innovations. It is based on an integral representation of the scattered field, a careful analysis on 
electromagnetic scattering and the behavior of the fundamental solutions. Numerically, it involves only 
computing inner products of the measured scattered field E^ with fundamental solutions to the Maxwell 
system located at sampling points over the measurement surface F. Hence it is strictly direct and does 
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not use any matrix operations or minimizations, and its implementation is very straightforward. Our 
extensive numerical experiments indicate that it can provide an accurate and reliable estimate of the 
scatterer support, even in the presence of a fairly large amount of noises in the data. Hence, it represents 
an effective yet simple computational tool for reliably detecting the scatterer support. In practice, a 
rough estimate of the scatterer support may be sufficient for many purposes [ ]. And if desired, one can 
obtain a more refined estimate of the medium scatterers by, e.g., any aforementioned indirect imaging 
methods or the one from [ ], using the estimate from the direct sampling method to provide an initial 
computational domain, which is usually much smaller than the originally selected sampling domain. Since 
indirect imaging methods often involve highly nonlinear optimization processes, an accurate and smaller 
initial domain can essentially reduce the entire computational efforts. 

The direct sampling method developed below uses a sampling strategy, and its flavor closely resembles 
MUSIC and the LSM (see [25, ] for overviews). However, it differs significantly from these two techniques. 
Firstly, it works with a few incident fields, whereas the latter two require the full map (multi-static 
response matrix/far-field operator) or data from sufficiently many incidents. Secondly, our method does 
not perform any matrix operations, e.g. eigendecomposition in MUSIC or solving ill-posed integral 
equations in the LSM. Hence, our method is computationally efficient. Lastly, the noise is treated directly 
via the inner product, which automatically filters out the noise contribution, and thus the method is highly 
tolerant to noises. 

The rest of the paper is organized as follows. In Section 2, we recall an integral reformulation of the 
Maxwell system, cf. (1), recently derived in [ ], which plays an essential role in the derivation of the direct 
sampling method. Then we develop the method in Section 3 in detail, where a preliminary analysis of its 
theoretical performance is also provided. In Section 4, we provide two- and three-dimensional numerical 
experiments to illustrate its distinct features, i.e., accuracy and robustness, for both exact and noisy data. 
Technical details of the implementation are provided in the appendix. 



2 Integral representation of Maxwell System 

In this part, we recall an equivalent formulation of the Maxwell system (1), which is fundamental to 
the derivation of the direct sampling method. We begin with the definition of the fundamental solution 
G{x,y) to the scalar Helmholtz equation, i.e., 

{-A-e)G{x,y)=S{x-y) 

where 5{x — y) is the Dirac delta function with the singularity located at G M^. We know that G(x, y) 
has the following representations (see, e.g., [ ]) 



G{x,y) = < 



-H','\k\x-y\), d = 2; 

1 ik\x-y\ 

d = 3. 



47r \x — y\^ 



where the function Hq^^ refers to Hankel's function of the first kind and zeroth order. Using the scalar 
function G{x^y) we can define a matrix-valued function ^{x^y) by 

^{x,y)=eG{x,y)I + D^G{x,y) (2) 

where / G R'^^^ is the identity matrix and denotes the Hessian of G. Then we can verify by some 
direct calculations that V • ^{x^y) = and 

V X V X ^(x, y) - k^^{x, y) = S{x - y)I, (3) 

where (and in the sequel) the actions of the operators V- and Vx on a matrix- valued function are always 
understood to be operated columnwise. Hence, the matrix ^{x^y) defined by (2) is a divergence-free 
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fundamental solution to the Maxwell system (1) in the homogeneous space R^. Using the fundamental 
solution ^{x^y)^ the total electric field E{x) can be represented by the following integral equation 

E{x) = E'^[ ^x,y){n^-l)E{y)dy. (4) 

We note that the fundamental solution ^{x^y) involves a non-integrable singularity at x = Hence care 
must be exerted when interpreting the integral, in the case that the point x lies within the domain Q [ ]. 
Next we let 7^ = n^ — 1, which precisely profiles the inhomogeneities of the media. In particular, the support 
of T] coincides with the scatterer support Cl. Furthermore, we introduce the function J = (n^ — 1)E^ that 
is the induced electrical current caused by the medium inhomogeneities. Then, the total electric field 
E{x) satisfies [ , Theorem 9.1] 

E{x) = E\x) ^k^ f G{x, y)J{y) dy^V^ f G{x, y)divyJ{y) dy. 

Upon noting the reciprocity relation \/xG{x^y) = —\/yG{x^y) and applying integration by parts to the 
second integral term on the right hand side, we arrive at the following equivalent integral equation 

E{x)- [ G{x,y)PJ{y)dy = E\x) 

jRd 

where the operator P is defined by 

P0 = k'^lcl) + grad(div0). 

Thus, by multiplying the equation with the coefficient r^, we obtain an integral equation for the generalized 
electric current J 

J{x)-r] I G{x,y)PJ{y)dy = r]E\x), xeVt. (5) 

The above equation was rigorously justified in suitable function spaces in [ ], and it is very convenient for 
solving inverse problems; see [ ] for the inverse source problem and [ ] for inverse medium scattering. 
Compared with the whole-space Maxwell system, the integral equation (5) is defined over the scatterer 
support Vt since the induced current J vanishes identically outside the scatterer support Vt = supp(77). 
This reduces greatly the computational domain, and hence brings significant computational convenience. 
We shall adopt the integral equation (5) for the forward scattering simulation, which can be discretized 
numerically by the mid-point quadrature rule (cf. Appendix A for details). 

3 Direct sampling method 

In this section, we develop a novel direct sampling method to determine the locations, the number and the 
shape of the scatterers/inhomogeneities in electromagnetic wave propagation. It is based on an analysis 
for electromagnetic scattering. Methodologically, it extends our earlier work on the acoustic scattering 
[21]. We shall also provide an analysis of the theoretical performance of the method by examining the 
behaviors of the fundamental solution. For the sake of convenience, we introduce the domain which 
is the domain enclosed by the circular /spherical measurement surface F, and use (•, •) for the real inner 
product on and the overbar for the complex conjugate. 

3.1 Direct sampling method 

In this part, we develop a novel direct sampling method (DSM). The derivation relies essentially on 
the following two basic facts. The first is the representation of the scattered electric field E^ using the 
fundamental solution <l> by (cf. (4)): 

E'{x)= [ ^{x,y)J{y)dy Vx G F. (6) 
Jn 
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The second fact is an important relation for the fundamental solution ^{x^y) (2) to the Maxwell 
system (1). For any two arbitrary sampling points Xp and Xq which lie inside domain but are away 
from its boundary F, we have the following approximation 

^(^(x, Xp)p, Xq)q)ds ^ k-^^{p, ^Xp, Xq)q) Vp, g G . (7) 

Next, we derive this crucial correlation. To do so, we first show the following lemma. 

Lemma 3.1. Let Xp and Xq be two distinct points lying inside the domain ftr- Then for any constant 
vectors p G and q gW^, there holds that 



(V X <l>(x, Xq)q X n, ^{x, Xp)p) - (V x <l>(x, Xp)p x n, ^{x, Xq)q)ds = -2i(p, ^{^{xp, Xq))q). (8) 



Proof. It is easy to verify directly that the identify (V x V x = V x V x ^{x^Xp)p holds for 

any constant vector p G M^. Using this identify and (3) we have 

V X V X (^{x,Xp)p - k^^{x,Xp)p = 8{x - Xp)p Vp G , (9) 

V xV x^{x,Xq)q-k^^{x,Xq)q = S{x - Xq)q VgGM"^. (10) 

Taking the conjugate of equation (10) yields 

V X V X <l>(x, Xq)q - /c^$(x, Xq)q = S{x - Xq)q. (11) 

By taking (real) inner products of equation (9) with ^{x^Xq)q and of equation (11) with ^{x^Xp)p^ then 
subtracting the two identities, we arrive at 

[{(^{x,Xq)q,V X V X (^{x,Xp)p) - {^{x,Xp)p,V x V x (^{x,Xq)q)]dx 

(12) 

= (P, ^{Xp, Xq)q) - (g, ^{Xp, Xq)p). 

Next we apply the following integration by parts formula 

[{^{x,Xq)q,V X V X ^{x,Xp)p) - {^{x,Xp)p,V x V x ^{x,Xq)q)]dx 

(V X <l>(x, Xq)q X n, ^{x, Xp)p) - (V x <l>(x, Xp)p x n, ^{x, Xq)q)ds 

to the left hand side of (12) to obtain 

(V X $(x, Xq)q X n, ^(x, Xp)p) — (Vx$(x, Xp)p x n, ^(x, Xq)q)ds 

= (P, ^(^p, Xq)q) - {q, ^{Xp, Xq)p). 
Now the real symmetry of the fundamental solution ^{x^y) leads directly to the identify (8). □ 



Next, note that on a circular curve/spherical surface F, we can approximate the left hand side of 
identity (8) by means of the Silver-Miiller radiation condition for the outgoing fundamental solution 
^{x^y) to the Maxwell system, i.e., 

V X <l>(x, Xp)p X n = ik^{x, Xp)p + h.o.t. 

Thus we have the following approximations: 

V X <l>(x,Xp)p X n ~ i/c<l>(x, Xp)j9, 

V X ^{x.,Xq)q X n ~ — i/c<l>(x, Xq)g', 
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which are vahd if the points Xp and Xq are not close to the boundary F. Consequently, we arrive at the 
following important approximate relation: 



/.{ 



{ik^{x,Xp)p,^{x,Xq)q) + {ik^{x,Xq)q,^{x,Xp)p)^ds ^ -2i{p,^{^{xp, Xq))q). 



Upon simplifying the relation, we arrived at the desired relation (7). 

The relation (7) leads us to a very important observation: the following integral over the measurement 
surface T 

(<!>(•, Xp)p, ^(•, Xq)q)L^r) = J (^(^, Xp)p, ^(x, Xq)q)ds (13) 

has a maximum if Xp = Xq and decays to as |xp — Xq| tends to oo, in view of the decay property of the 
fundamental solution ^{x^y). 

These two basic facts (6) and (7) are crucial to the derivation of our direct sampling method. Next 
we consider a sampling domain 17 C l^r that contains the scatterer support 17, and divide it into a set of 
small elements {tj}. Then by using the rectangular quadrature rule, we arrive at the following discrete 
sum representation 

E^ix) = f <i>ix,y)Jiy)dy « ^ y,) J(y,) |r,|, (14) 
Jn 

where the point yj lies within the jth element r^, and \rj\ denotes the volume/area of the element tj. 
Noting the fact that the generalized induced current J vanishes identically outside the support 17, the 
summation in (14) is actually only over those elements that intersect with Q. We point out that, in 
practice, the scatterer support Q may consist of several disjoint subregions, each being occupied by a 
(possibly different) physical medium. By elliptic regularity theory [11, 31], the induced current J = r] E 
is smooth in each subregion, and thus according to classical approximation theory, the approximation in 

(14) can be made arbitrarily accurate by refining the elements {tj}. Nonetheless, we reiterate that the 
relation (14) serves only the goal of motivating our method, and it is not needed in the implementation 
of our method. Physically, the relation (14) can be interpreted as follows: the scattered electric field 

at any fixed point x G F is a weighted average of that due to the point scatterers located at {yj} lying 
within the true scatterer Q. _ 

The approximate relations (7) and (14) together give that for any sampling point Xp e Q and any 
constant vector G M^, there holds 

{E',^-,Xp)q)L^r) - {^^{■,yj)J{yj)\Tjl^{-,Xp)q\ 

\ ^ I L2(r) 

3 

«fc-i^|r,|(J(y,),9($(xp,y,))q). 

3 

Due to the behavior of the fundamental solution (see the discussions in Section 3.2), the approximation 

(15) indicates that the inner product {E^ ^^{•^Xp)q) may take significant values if the sampling point Xp 
is close to some physical point scatterer yj and small values otherwise, thereby acting as an indicator 
function to the presence of scatterers and equivalently providing an estimate of the scatterer support Vt. 
In summary, these observations lead us to the following index function: 

^fixj)]q) = r— — n-^ , Mx^eil. (16) 

^ ^'^^ \\E^\\L2ir)m-.Xp)q\\L^^r) ' ^ ^ 

Here ^{■^Xp)q acts as a probe/detector for the scatterers. In principle, the choice of the polarization 
vector q in the index ^{xp]q) can be quite arbitrary. Naturally, it is expected that the choice of q will 
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affect the capability of the function ^{x^Xp)q for probing scatterers. The analysis in Section 3.2 will 
shed insights into the probing mechanism and will provide useful guidelines for the choice of q. Our 
experimental experiences indicate that the choice q = i.e., the polarization of the incident field 
works very well in practice. We note that, apart from providing an estimate of the scatterer support, the 
index ^(x^; q) provides a likelihood distribution of the inhomogeneities in 17, which up to a multiplicative 
constant may be used as the an initial estimate of the distributed coefficient r] [ ]. 

The DSM is one sampling type imaging technique. There are several existing sampling techniques, e.g., 
MUSIC and the LSM. The DSM differs from these techniques in several aspects. First, the DSM involves 
computing only inner products of the measured scattered field with the probing functions ^{x,Xp)q, 
i.e., the fundamental solutions with singularity concentrated at sampling points. Therefore, it does not 
involve any expensive linear algebraic operations, like solving first-kind ill-posed integral equations in the 
LSM and eigenvalue decomposition for computing the signal space (noise space) in the MUSIC. Secondly, 
the DSM is applicable to a few (e.g., one or two) incident fields, whereas the LSM requires the far 
field map and the MUSIC requires the multi-statistic matrix, which necessitates multiple measurements. 
Therefore, the DSM is particularly attractive in the case of a few scattered data. Thirdly, the DSM treats 
noise directly via the inner product. The noise roughly distributes the energy in all frequencies equally, 
whereas the fundamental solutions are smooth with their energy concentrated on low frequency modes. 
Therefore, the high-frequency modes in the noise are roughly orthogonal to the (smooth) fundamental 
solutions. In other words, the presence of noise does not cause much changes in the inner product, but 
the normalization. As a consequence, the DSM is highly tolerant to data noise. 

3.2 Analysis of the index function 

In this part, we analyze the theoretical performance and the mechanism of the DSM by analytically 
and numerically studying the fundamental solution ^{x^y) in (2) and the index ^ in (16) for one single 
point scatterer. First, we recall the crucial role of the approximate relation (7): the fundamental solution 
^(xp, Xq) is nearly singular, and assumes very large values for Xp close to Xq. More precisely, the extremal 
property of ^{^{xp^Xq)) paves the way for accurate support detection. This observation is evident for 
the scalar Helmholtz equation in view of the identity ^{G{xp^Xq)) — \ J(^(k\xp — Xq|) (in two-dimension), 
where Jo is the Bessel function of the first kind of order zero. However, for the Maxwell system, the 
fundamental solution ^{x^y) contains 4 or 9 different entries, and each term exhibits drastically different 
behavior. So with each different polarization vector q^ the probing function <I>(x, Xp)q^ as well as the index 
function ^{xp\q) in (16), mixes different components together and may yield quite different behaviors. 
Next we shall investigate more closely the properties of ^{x^y) in order to have a better understanding 
and interpretation of the index function ^[xp\q) in (16). 

We begin with an important observation on the trace of the fundamental solution ^{x^y). 

Proposition 3.1. For d = 2,3^ the fundamental solution ^{x^y) satisfies the following relation 



Proof. Let us begin with the two-dimensional case, i.e., transverse electric mode, i.e., E = (£^i,£^2,0) 
and H = (0,0, i^a). Then the Maxwell system (1) should be understood as 



tT^x,y) = {d-l)eG{x,y). 



(17) 



icjfio = (V X E)s 



dE2 
dxi 



dEi 

dX2 




The fundamental solution G{x^y) is given by G{x^y) = ^HQ^\k\x — y\). To evaluate the Hessian 
D'^G{x^y), we first recall the recursive relation for the derivative of Hankel functions [±] 
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Hence, by chain rule and product rule, we deduce that 



H},'\k\x-y\)=eH!,'>{k\x-y\) 



dxidxj 



^ {x-y)i{x-y)j 
\x-y\'^ 



fc^J?W(fc|a:-y|)5i,,-, 



where 5ij is the Kronecker delta function. Consequently, the fundamental solution ^{x^y) is given by 
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H^'\k\x-y\) 



H['\k\x-y\) 
k\x-y\ 



{x-y)i{x-y)j 
\x-y\'^ 



Upon noting the recursive relation ^^^^ = H^\z) + H2'\z) for Hankel functions [ ], we deduce that 



the trace of the fundamental solution satisfies 



ii{^{x,y)) 



2H^^\k\x-y\)-2 



H[^\k\x-y\) 
k\x-y\ 



■H^^\k\x-y\) 



'-^Hi^\k\x-y\) = eG{x,y). 



Next we turn to the three-dimensional case, i.e., G{x^y) 



\k\x — y\ 



J_e 

47r \x-y\ 



A direct calculation yields 



^{x,y)i 



--G{k\x-] 



jx-y) 
I 



^x-y)i{x-y) 
\x-y\^ 



y)i{x-y)j \ _ fc _ 3{x-y)i{x-y)j \ 1 
x-y\^ J y^^d \x-y\^ J \x-y\^ 

J \x-y\_ 



Consequently, tr(<l>(x, ?/)) = 2k'^G{x^y). 



□ 



It follows from (17) that the sum of the diagonal components of the fundamental solution ^{x^y) 
to the Maxwell system behaves like the scalar fundamental solution G{x^y)^ which, according to the 
experimental evidences in [21], can provide an excellent probing function of the scatterer support. To be 
more precise, this yields 



k~^^{tT^{x,y)) 



k^Jo{k\x-y\) 



2k 



2 sin(/c|x 



(18) 



k\x - y\ 



d = 3. 



Therefore, the crucial quantity /c~-^^(tr<l>(xp, Xq)) attains its maximum at Xp = x^, indicating the presence 
of the scatterers. We note that apart from the global maximum, there are also a number of local maxima, 
which introduce "ripples" into the resulting index. However, we note that the positions of the local 
maxima depend only on the known wavenumber /c, and can be strictly calibrated. 

Next, we examine combinations of the components of the fundamental solution ^{x^y). We first note 



that each independent component ^^ij{x^Xq) does have a significant maximum at x 



like (18). 



However, the magnitude of each component differs dramatically. In Fig. 1, we show the inner products 
of the components <l>ii, <l>22, ^12 with themselves (termed cross product for short) and the diagonal sum 
as a function of the sampling point Xp G 1^, i.e., |(<l>ii(xp, x), <I>ii(xg, x)/,2,^r)) | etc. We note that these 
cross products form the building blocks of the index function ^(xp;g), and their behaviors completely 
determine the performance of the DSM. The cross products for <I>ii and ^22 exhibit significant directional 
resonance effects but their sum, i.e., (<I>ii(x, Xg), <^>ii(x, Xp))l2^y) + (^22 (x, ^q), *^*22(x, Xp))l'^{t) does only 



assume one significant maximum at Xq 



without much resonance effects. In practice, due to the 



mixing of these different components in the index ^{xp]q)^ the performance of the DSM may not be as 
good as that for the scalar Helmholtz case in the case of the scattered electric field of one polarization 
p. Nonetheless, with multiple polarizations pk and all components of the respective scattered field ^ 
one might be able to approximate the trace tr<I>(x, ?/), and (18) can be applied. If this is indeed the case, 
then the performance of the sampling method would equal that for the scalar case. 
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|(<3^i2(x,Xp),<l>i2(x,Xg))i^2(r)| diagonal sum 



Figure 1: Cross products of <l>ii, ^22^ ^12, with the point scatterer located at the Xq = (— 1,0), and the 
diagonal sum, over the sampling domain Q = [—2,2]^. The quantities are normalized. 

Index function for multiple incident fields. To arrive at an applicable effective index function for 
multiple incidents, we examine the behavior of the two-dimensional fundamental solution ^{x^y) more 
closely. We first consider the incident direction di = ^(1, 1) and the polarization pi = ^(1,-1), and 
assume that the local induced electrical current J is proportional to the polarized incident wave E'^ within 
the inhomogeneity located at Xq. Upon ignoring the normalization, the index function "^{xp'^q) with the 
choice q = Pi is roughly determined by 

{^{x,Xq)pi,^{x,Xp)pi)L2^r) ={^ll{x,Xq) - ^i2{x, Xq) , ^ii{x, Xp) - <^>12 ^p) ) L2 (F) 

+ {^12{x,Xq) - ^22{x,Xq),^i2{x,Xp) - ^22{x , Xp)) l^F) , 

where each term in the inner products, like ^ii{x^Xp) — ^i2{x^Xq) and so on, corresponds to one com- 
ponent of the vector field E^. Similarly, for the incident direction d2 = ^(1,-1) and the polarization 
P2 = ^(1,1), the index function ^{xp; q) with the choice q = P2 is roughly determined by 

{^{x,Xq)p2,^{x,Xp)p2)L^r) ={^ll{x,Xq) + ^I2{x, Xq) , ^ii{x, Xp) + <I>12 (x, Xp) ) ^^2 (r) 

+ {^12{x,Xq) + ^22{x,Xq),^i2{x,Xp) + ^22 ^p))L2(r) • 

The quantities \{^{x^Xq)pi^^{x^Xp)pi)L2(^Y)\ and \{^{x^Xq)p2^^{x^Xp)p2)L^(^r)\ after normalization are 
shown in Fig. 2: they exhibit strong resonances. This is attributed to the ripples predicted from the 
relation (18) and mixing of different components of the fundamental solution <l>. Note that the resonances 
show up at different regions, with their locations depending on the incident direction/polarization. We 
reiterate that in case of one single point scatterer, the resonance is completely predictable and one 
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might be able to remove them from the reconstructed images with suitable postprocessing. However, 
in the presence of multiple scatterers, the discrimination of the resonances from the true scatterers can 
be highly nontrivial. Therefore, it is highly desirable to design a sampling method that is free from 
significant resonances ab initio, by suitably adapting the choice of q in the probing function ^{x^Xp)q 
and combining the indices for multiple polarizations. Let us consider the two polarizations pi = ^(1,-1) 
and p2 = ^(1, 1), and measure all components of the scattered electric field E^, in the hope of focusing 
on the genuine inhomogeneities. Then the performance of the combined index with pi = ^(1, —1) and 
P2 = ^(1, 1) is approximately given by 

{^{X, Xq)pi,^{x, Xp)pi)L^r) + (^(^, Xq)p2,^{x, Xp)p2) L^F) 
= (<^>ll(x, Xq), <^>il(x, Xp))i^2(r) + 2{^12{X, Xq), ^I2{x, Xp))i^2(r) + {^22{X, Xq), <^>22(^, ^p))L2(r) • 

Even though the indices for the polarizations pi and p2 separately both have their own shadows of the 
point scatterer at Xq = ( — 1,0) (see Fig. 2(a)-(b)), the combination can isolate the point scatterer very 
well (see Fig. 2(c)). Hence, the combination of multiple polarizations can remedy the undesirable rippling 
phenomena. This suggests us the following combined index ^c(^p) for the multiple polarizations {pi}f'=i : 

1 ^ 

^c(%) = ^^^(^p;pO 

where "^{xp'^pi) refers to the index function for the polarization pi (hence the Ith incident field). 




(a) \{^{x,Xp)pi,^{x,Xq)pi)L^r)\ (b) \{^{x, Xp)p2,^{x, Xq)p2) L^r)\ (c) combination 

Figure 2: The cros^s products for polarizations pi, p2 and the combination with direct sum over the 
sampling domain Q = [—2,2]^. The point scatterer is located at Xq = (— 1,0). The quantities are 
normalized. 



4 Numerical experiments and discussions 

In this section, we present two- and three-dimensional numerical examples to showcase the features, e.g., 
the accuracy and robustness, of the proposed DSM for detecting scatterers. The examples are designed 
to illustrate the features of the method and to validate the theoretical findings in Section 3. Numerical 
results for both exact and noisy data will be presented. In all examples, the wavelength A is set to 1, and 
the wavenumber k is 27r. The noisy scattered electric data Eg is generated pointwise by 

El{x) = E'{x)^emdix\E'{x)\C{x), 

where e denotes the relative noise level, and both real and imaginary parts of the random variable ({x) 
follow the standard Gaussian distribution with zero mean and unit variance. All the computations were 
performed on a dual-core laptop using MATLAB R2009a. 
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4.1 Two-dimensional examples 

Now we present numerical results for two-dimensional examples. Unless otherwise specified, two incident 
fields, i.e., di = ^(1,1)* and d2 = ^(—1,1)* (accordingly, the polarizations pi = ^(1,-1)* and 
P2 = ^(1,1)*), are employed. For each incident field the scattered electric field is measured at 
30 points uniformly distributed on a circle of radius 5 centered at the origin. The sampling domain Q is 
fixed at [—2, 2]^, which is divided into small squares of equal side length h = 0.01. The index function ^ 
as an estimate to the scatterer support will be displayed for each example. 

Our first example shows the method for one single scatterer, which confirms the theoretical analysis 
of the index in Section 3.2. 

Example 1. The example considers one square scatterer of side length 0.3 centered at the point (— 1,0). 
The true inhomogeneity coefficient r]{x) of the scatterer is set to be 1, namely n{x) = a/2. 




(a) true scatterer (b) index (c) index ^(xp;p2) (d) index 



Figure 3: Numerical results for Example 1. The first and second row refers to index functions for the 
exact data and the noisy data with e = 20% noise. 

The numerical results for Example 1 are shown in Fig. 3. We observe that for sampling points close to 
the physical scatterer, the index function ^ is relatively large; otherwise it takes relatively small values. 
Note that the image for one incident field E'^ exhibits obvious resonances, with resonance locations 
depending on the incident direction. Here, the resonance behavior agree excellently with the theoretical 
analysis for one single point scatterer in Section 3.2, hence it might be removed by applying a suitable 
postprocessing procedure. Nonetheless, the use of two incident fields can greatly mitigate the resonances. 
Therefore, the index function does provide an accurate and reliable indicator for the location of the 
scatterer; see Fig. 3(d). The presence of e = 20% noise in the measured data does not affect the shape of 
the index function thereby showing the robustness of the method with respect to noise. 

Our second example illustrates the method for two separate scatterers. 

Example 2. We consider two square scatterers, with the inhomogeneity coefficient r] being 1 in both 
scatterers. The following two cases are investigated: 

(a) The two scatterers are of side length 0.2, and located at (—0.8,-0.7) and (0.3,0.8); respectively. 

(b) The two scatterers are of side length 0.3^ and located at (—0.45, —0.35) and (0.05,0.15); respectively. 
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(a) true scatterers (b) index (c) index '^{xp;p2) (d) index 

Figure 4: Numerical results for Example 2(a). The first and second row refer to index functions for the 
exact data and the noisy data with e = 20% noise. 



The two scatterers in Example 2(a) are well apart from each other. In this example, each of the two 
incident fields tends to highlight only one of the two scatterers, with the index value for one scatterer 
being much higher than for the other; see Figs. 4(b)-(c). Since the two scatterers are well apart, the 
interactions between them is somehow weak, and the resonance pattern for the point scatterer is well 
kept. However, the two incident fields together give a clear discrimination of the two scatterers, with 
their locations satisfactorily recovered for both exact data and the data with e = 20% noise; see Fig. 4(d). 




(a) true scatterers 



(b) index ^(xp;pi) 



(c) index ^{xp]p2) 



(d) index 



Figure 5: Numerical results for Example 2(b). The first and second row refer to index functions for the 
exact data and the noisy data with e = 20% noise. 



The two scatterers in Example 2(b) stay very close to each other. We observe that, for the incident 
direction 6/2, apart from the strong resonances, the locations for the scatterers cannot be directly inferred 
since the index function ^ relates the two scatterers into an elongated ellipse shape. Nonetheless, the 
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resonances were almost completely removed from the estimate when using two incident fields. Conse- 
quently, the estimate of the locations of the scatterers is very impressive: the two scatterers are still well 
separated despite their closeness, with their locations correctly estimated, for up to e = 20% noise in the 
data. Although not presented, we would like to remark that in the case of very high noise levels, e.g., 
6 = 40%, the estimate tends to connect the two scatterers, and also some spurious modes have emerged. 
Our next example considers the more challenging case of three neighboring scatterers. 

Example 3. This example consists of three neighboring square scatterers of width 0.15.* one located at 
(— |, — |)^ one located at {—^^—^), one located at (— |). The inhomogeneity coefficients of all 
three scatterers are set to be r] = 1. 




(a) true scatterers (b) index (c) index '^{xp;p2) (d) index 



Figure 6: Numerical results for Example 3. The first and second row refer to index functions for the 
exact data and the noisy data with e = 20% noise. 

In this example, the three scatterers stay very close to each other, especially the upper two, thus it 
is numerically very challenging to separate them. This is also reflected in the fact that each individual 
incident fleld tends to combine two of the three scatterers into one larger chunk, which is true for both 
the exact data and noisy data; see Figs. 6 (b)-(c). Therefore, it is diflicult to tell from either Fig. 6(b) 
or Fig. 6(c) the number and locations of the scatterers. The latter is effectively remedied by using two 
incident fields together; see Fig. 6, where the three scatterers are vividly separately from each other. 
However, the interactions between the scatterers focus the strength on the scatterer to the right, and 
diminish slightly the strength of the scatterer to the left. 

Next we consider a ring-shaped scatterer. 

Example 4. This scatterer is one ring-shaped square scatterer located at the origin, with the outer and 
inner side lengths being 0.6 and OA, respectively. The coefficient r] of the scatterer is 1. 

The ring-shaped scatterer represents one of the most challenging objects to resolve, and it is highly 
nontrivial even with multiple data sets. The results with the exact data and e = 20% noise in the data 
are shown in Fig. 7. It is observed that with just two incident waves, the method can provide a quite 
reasonable estimate of the ring shape, and it remains very stable for up to e = 20% noise in the data. 

4.2 Three-dimensional example 

The last example shows the feasibility of the method for three-dimensional problems. 
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(a) true scatterers (b) index (c) index '^{xp;p2) (d) index 



Figure 7: Numerical results for Example 4. The first and second row refers to index functions for the 
exact data and the noisy data with e = 20% noise. 

Example 5. We consider two cubic scatterers with side length 0.2; one centered at (0.4,0.3,0.3) and the 
other at (—0.4,0.3,0.3); respectively, and the coefficient 77 in both scatterers is taken to be 1. 

For this example, we take two incident fields, with the incident directions di = d2 = 1? 1)* 

the polarization vectors pi = ^(1, —2, 1)* and p2 = ^(1^ 1? —2)*. The scattered field is measured at 
the points on a uniformly distributed mesh of 10 x 10 on each face of the cube of edge length 10. The 
sampling domain Q for evaluating the index functions is taken to be [—2, 2]^. The problem is discretized 
with a mesh size 0.02. The numerical results are shown in Fig. 8. We observe that the support estimated 
by the index ^ agrees very well with the exact one, the magnitude of the index ^ decreases quickly away 
from the boundary of the true scatterers. The presence of e = 20% data noise (cf. Fig. 9), seems to cause 
no obvious deterioration of the accuracy of the index ^ as compared with the exact data. By examining 
the cross-sectional images of the index, we found that in comparison with the two-dimensional problems, 
the rippling phenomenon seems less pronounced for this three-dimensional example. 

5 Conclusions 

We have developed a novel direct sampling method for the inverse electromagnetic medium scattering 
problem of estimating the support of inhomogeneities from scattered electric near-field data. It was 
derived based on an integral representation of the scattered field via the fundamental solution, a careful 
analysis on electromagnetic scattering and the behavior of the fundamental solutions. The method is 
particularly attractive in that it is applicable to a few incident fields. Methodologically, it involves only 
computing inner products of the scatted electric field with fundamental solutions located at the sampling 
points, hence it is strictly direct, straightforward to implement, computationally very efficient, and very 
robust with respect to the presence of data noise. The experimental results for two- and three-dimensional 
examples indicate that it can provide, after thresholding with a suitable cutoff value, a quite satisfactory 
estimate of the shapes of the scatterers from the measured data corresponding to only one or two incident 
directions, even in the presence of a large amount of noise. 

In the present study, we have focused only on the very cheap direct sampling method for the re- 
constructions of the scatterers. It is natural to enhance the reconstructions by other more expensive 
but more accurate methods. More specifically, one can use the reconstructions by the direct sampling 
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(c) index ^(xp;p2) (d) index 

Figure 8: Numerical results for Example 5 with exact data. 



method as the initial computational domain for those more refined methods to retrieve more accurate 
shapes of the scatterers and their physical medium properties, e.g., Tikhonov regularization [ ]. Since 
the indirect imaging methods often involve highly nonlinear optimization processes, a more accurate and 
much smaller initial domain can essentially reduce the relevant computational costs. 

It is also interesting to see the extensions of the direct sampling method for other important scattering 
scenarios, e.g., scattering from lines (cracks), far-field measurements and multiple-frequency data, as well 
as their theoretical justifications. 
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A Numerical method for forward scattering 

In this part we describe our numerical method for the integral equation (5) for two-dimensional domains 
ft. We denote by J the index set of grid points Xj = (^j^, ^^3)' ^ ~ Oi'^s) G JJ, of a uniformly distributed 
mesh with a mesh size h > and the square cells Bj given by 

~ ^jid2 ~ ^^311^32) ^ [~ 2 ' 2"] ^ [~2' 2"] 
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(c) index ^(xp;p2) (d) index 

Figure 9: Numerical results for Example 5 with e = 20% data noise. 



for every tuple j = (ji,j2) in the index set JJ. Further, we assume that the set Uj^^Bj contains the 
scatterer support ^. We shall approximate the integral operator in (5) by the mid-point quadrature rule, 
i.e., 

Jk^Vk Y,Gk,j{PJ)jh^ =VkE\xk) (19) 

jei 

where Jk = J{xk) and rjk = Tjixk)^ and the off-diagonal entries Gkj and the diagonal entries Gk^k are 
given by G^j = G{xk,Xj) and 

Gk,k = T2 I G{x,0)dx, 

respectively. The diagonal entries Gk,k can be accurately computed by a tensor-product Gaussian quadra- 
ture rule. To arrive at a fully discrete scheme, we further approximate the crucial term PJ in equation 
(19) by the central finite difference scheme: 

where Dx^xj refers to taking the second-order derivative with respect to Xi and Xj. In practice, we shall 
approximate Dx^xj by central finite difference scheme, i.e., 

Dx,x^ =H^I, Dx^x2 = I®H, Dx,x2 =D^D, 

where (g) is the Kronecker product for matrices, H and D are the tridiagonal matrices for the second 
derivative and the first derivative, respectively. The resulting system can be solved using standard 
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numerical solvers, e.g., Gaussian elimination, if the cardinality of the index set JJ is medium, or iterative 
solvers like generalized minimal residual method (GMRES) [ ] should be applied if the cardinality of 
JJ is large. Clearly, the extension of the procedure to 3D problems is straightforward. Our numerical 
experiences indicate that tens of GMRES iterations already yield a very accurate solution to the linear 
system. 
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